En cada sección podrás ejecutar y analizar los comandos básicos para leer, visualizar y realizar los análisis básicos utilizados en Dendrocronología. La mayoría de análisis que se llevaran acabo son iguales o muy similares a los que se realizan en los programas originales programados en FORTRAN llamados COFECHA y ARSTAN (https://www.ldeo.columbia.edu/tree-ring-laboratory/resources/software).
En este caso utilizaremos el ambiente de R y utilizaremos el paquete dplR programado y desarrollado por Andy Bunn y Mikko Korpela en el 2014 pueden bajar el articulo dándole clic al siguiente enlace. https://drive.google.com/uc?export=download&id=1ps7tZ0RTwuIr2P_h-ugHkbgAs3Jrdxaa
La idea y objetivo principal de estos ejercicios es que puedas familiarizarte con este tipo de ambiente, y puedas aprender haciendo. A lo largo de este documento hay un set de comandos y ejercicios que podrás ejecutar y modificar con tus propios archivos y o mediciones. Si lo modificaste y necesitas el archivo original lo puedes clonar de nuevo desde la pagina de github. https://github.com/paulszejner/Intro_dendrocrono.git
A continuación vamos a instalar y activar el dplR y se conectaran a la base de datos disponibles del international Tree Ring Data Bank (ITRDB) mediante la web
El primer paso de buenas costumbres y buenas practicas es iniciar con el ambiente de R limpio utilizando el comando rm(list=ls()) podrás limpiar todos los objetos que pueden estar activos en tu consola.
# para instalar el paquete de inertes hay que borrar el numeral que aparece antes del "install.packages(dplR)" y ejecutar esa linea
# install.packages(dplR)
library(dplR) # librería para análisis dendrocronológicos
Dat_ringwidth <- read.rwl("https://www.ncei.noaa.gov/pub/data/paleo/treering/measurements/northamerica/mexico/mexi042.rwl")
Attempting to automatically detect format.
Assuming a Tucson format file.
There appears to be a header in the rwl file
There are 20 series
1 RAN01A 1950 2006 0.001
2 RAN01B 1929 2006 0.001
3 RAN02A 1957 2006 0.001
4 RAN02B 1933 2006 0.001
5 RAN03A 1929 2006 0.001
6 RAN03B 1942 2006 0.001
7 RAN04A 1921 2006 0.001
8 RAN05A 1925 2006 0.001
9 RAN05B 1934 2006 0.001
10 RAN06A 1925 2006 0.001
11 RAN07A 1924 2006 0.001
12 RAN07B 1918 2006 0.001
13 RAN08A 1913 2006 0.001
14 RAN08B 1913 2006 0.001
15 RAN09A 1913 2006 0.001
16 RAN09B 1913 2006 0.001
17 RAN10A 1865 2006 0.001
18 RAN11B 1946 2006 0.001
19 RAN12A 1936 2006 0.001
20 RAN12B 1936 2006 0.001
Dat_Earlywood <- read.rwl("https://www.ncei.noaa.gov/pub/data/paleo/treering/measurements/northamerica/mexico/mexi042e.rwl")
Attempting to automatically detect format.
Assuming a Tucson format file.
There appears to be a header in the rwl file
There are 20 series
1 RAN01A 1950 2006 0.001
2 RAN01B 1929 2006 0.001
3 RAN02A 1957 2006 0.001
4 RAN02B 1933 2006 0.001
5 RAN03A 1929 2006 0.001
6 RAN03B 1942 2006 0.001
7 RAN04A 1921 2006 0.001
8 RAN05A 1925 2006 0.001
9 RAN05B 1934 2006 0.001
10 RAN06A 1925 2006 0.001
11 RAN07A 1924 2006 0.001
12 RAN07B 1918 2006 0.001
13 RAN08A 1913 2006 0.001
14 RAN08B 1913 2006 0.001
15 RAN09A 1913 2006 0.001
16 RAN09B 1913 2006 0.001
17 RAN10A 1865 2006 0.001
18 RAN11B 1946 2006 0.001
19 RAN12A 1936 2006 0.001
20 RAN12B 1936 2006 0.001
Dat_latewood <- read.rwl("https://www.ncei.noaa.gov/pub/data/paleo/treering/measurements/northamerica/mexico/mexi042l.rwl")
Attempting to automatically detect format.
Assuming a Tucson format file.
There appears to be a header in the rwl file
There are 17 series
1 RAN02A 1960 2006 0.001
2 RAN02B 1933 2006 0.001
3 RAN03A 1929 2006 0.001
4 RAN03B 1942 2006 0.001
5 RAN04A 1921 2006 0.001
6 RAN05A 1925 2006 0.001
7 RAN05B 1934 2006 0.001
8 RAN06A 1925 2006 0.001
9 RAN07A 1924 2006 0.001
10 RAN07B 1918 2006 0.001
11 RAN08A 1913 2006 0.001
12 RAN09A 1913 2006 0.001
13 RAN09B 1913 2006 0.001
14 RAN10A 1865 2006 0.001
15 RAN11B 1946 2006 0.001
16 RAN12A 1936 2006 0.001
17 RAN12B 1936 2006 0.001
# También se puede bajar desde los FPTs del NOAA
# Dat_ftp <- read.rwl("ftp://ftp.ncdc.noaa.gov/pub/Data/paleo/treering/measurements/northamerica/mexico/mexi042.rwl")
str(Dat_ringwidth) # la función str te da la estructura del objeto llamado Datos
Classes ‘rwl’ and 'data.frame': 142 obs. of 20 variables:
$ RAN01A: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN01B: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN02A: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN02B: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN03A: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN03B: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN04A: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN05A: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN05B: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN06A: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN07A: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN07B: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN08A: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN08B: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN09A: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN09B: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN10A: num 1.08 1.47 1.7 1.08 1.48 ...
$ RAN11B: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN12A: num NA NA NA NA NA NA NA NA NA NA ...
$ RAN12B: num NA NA NA NA NA NA NA NA NA NA ...
### Grafica los datos para el control de calidad -------------------------------
plot.ts(Dat_ringwidth[,20])
plot.ts((Dat_ringwidth[,2]), col="darkred")
plot.ts((Dat_ringwidth[,3]), col="darkred")
plot.ts((Dat_ringwidth[,6]), col="darkred")
Podemos revisar todos los detalles del archivo con el cual vamos a trabajar
# cuantas series de tiempo están en este archivo? preguntando cuantas columnas tenemos es una forma de responder la pregunta
n_series <- ncol(Dat_ringwidth)
# ahora hacemos un plot con todas las serias en su estado crudo o RAW sin ningún tratamiento alguno
plot.ts(Dat_ringwidth[,20])
for(s in 1:n_series){lines(Dat_ringwidth[,s]) } # esta forma de plotear todas las series es un for loop usando cada columna en cada uno de los loops
dplR llamado rwl.report()rwl.report(Dat_ringwidth)
Number of dated series: 20
Number of measurements: 1618
Avg series length: 80.9
Range: 142
Span: 1865 - 2006
Mean (Std dev) series intercorrelation: 0.4451521 (0.08958034)
Mean (Std dev) AR1: 0.7406 (0.08279931)
-------------
Years with absent rings listed by series
None
-------------
Years with internal NA values listed by series
None
En dplR hay otras opciones para graficar los datos que tenemos. Por ejemplo con mas orden y con algunas características que pueden ser útiles. La funcion que se utiliza es el famoso spaguetti plot spag.plot funcion del dplR
spag.plot(Dat_ringwidth, zfac = 1)
spag.plot(Dat_ringwidth, zfac = 3)
spag.plot(Dat_ringwidth, zfac = 0.3)
spag.plot(Dat_Earlywood, zfac = 0.3)
spag.plot(Dat_latewood, zfac = 0.3)
Luego de fechar visualmente las muestras de madera, se procede a su medición, Durante la medición se pueden cometer algunos errores, lo cual puede perjudicar la señal común de los árboles. Por lo tanto debemos corroborar las correlaciones entre todas las series antes de proceder a cualquier análisis.
#simple correlación entre series y series con su periodo en común
cor_matrix <- round(cor(Dat_ringwidth, use = "pairwise"), digits = 2) # la función cor va a utilizar una matriz o un data frame para calcular las correlaciones entre todas las columnas presentes en la base de datos. generando una matriz cuadrada con le numero de columnas y filas igual al número de series.
cor_matrix
RAN01A RAN01B RAN02A RAN02B RAN03A RAN03B RAN04A RAN05A RAN05B RAN06A RAN07A RAN07B RAN08A RAN08B RAN09A RAN09B RAN10A
RAN01A 1.00 0.54 0.19 0.50 0.56 0.50 0.35 0.70 0.78 0.67 0.47 0.42 0.79 0.80 0.60 0.56 0.74
RAN01B 0.54 1.00 0.16 0.24 0.07 0.05 0.24 0.34 0.30 0.28 -0.18 -0.20 0.51 0.51 0.09 0.07 0.35
RAN02A 0.19 0.16 1.00 0.52 0.16 0.05 0.28 0.26 0.30 0.11 0.27 0.15 0.22 0.28 0.17 0.15 0.24
RAN02B 0.50 0.24 0.52 1.00 0.32 0.41 0.22 0.37 0.45 0.53 0.35 0.32 0.33 0.40 0.55 0.44 0.35
RAN03A 0.56 0.07 0.16 0.32 1.00 0.79 0.41 0.53 0.68 0.51 0.61 0.56 0.52 0.58 0.55 0.62 0.58
RAN03B 0.50 0.05 0.05 0.41 0.79 1.00 0.47 0.58 0.73 0.80 0.58 0.62 0.55 0.61 0.66 0.79 0.61
RAN04A 0.35 0.24 0.28 0.22 0.41 0.47 1.00 0.63 0.48 0.56 0.18 0.34 0.37 0.51 0.58 0.61 0.34
RAN05A 0.70 0.34 0.26 0.37 0.53 0.58 0.63 1.00 0.72 0.63 0.06 0.23 0.70 0.76 0.50 0.59 0.71
RAN05B 0.78 0.30 0.30 0.45 0.68 0.73 0.48 0.72 1.00 0.83 0.41 0.42 0.75 0.82 0.60 0.63 0.79
RAN06A 0.67 0.28 0.11 0.53 0.51 0.80 0.56 0.63 0.83 1.00 0.20 0.43 0.59 0.68 0.67 0.69 0.60
RAN07A 0.47 -0.18 0.27 0.35 0.61 0.58 0.18 0.06 0.41 0.20 1.00 0.66 -0.04 0.05 0.45 0.38 0.10
RAN07B 0.42 -0.20 0.15 0.32 0.56 0.62 0.34 0.23 0.42 0.43 0.66 1.00 0.25 0.28 0.62 0.65 0.31
RAN08A 0.79 0.51 0.22 0.33 0.52 0.55 0.37 0.70 0.75 0.59 -0.04 0.25 1.00 0.93 0.35 0.46 0.66
RAN08B 0.80 0.51 0.28 0.40 0.58 0.61 0.51 0.76 0.82 0.68 0.05 0.28 0.93 1.00 0.45 0.51 0.60
RAN09A 0.60 0.09 0.17 0.55 0.55 0.66 0.58 0.50 0.60 0.67 0.45 0.62 0.35 0.45 1.00 0.89 0.35
RAN09B 0.56 0.07 0.15 0.44 0.62 0.79 0.61 0.59 0.63 0.69 0.38 0.65 0.46 0.51 0.89 1.00 0.51
RAN10A 0.74 0.35 0.24 0.35 0.58 0.61 0.34 0.71 0.79 0.60 0.10 0.31 0.66 0.60 0.35 0.51 1.00
RAN11B 0.33 0.24 0.44 0.31 0.01 -0.14 -0.32 0.13 0.14 -0.10 -0.02 -0.15 0.27 0.17 -0.27 -0.30 0.25
RAN12A 0.41 0.15 0.28 0.29 0.37 0.58 0.25 0.32 0.48 0.48 0.08 0.09 0.39 0.48 0.24 0.29 0.36
RAN12B 0.63 0.33 0.20 0.24 0.45 0.53 0.17 0.51 0.67 0.47 0.02 0.12 0.70 0.70 0.17 0.24 0.67
RAN11B RAN12A RAN12B
RAN01A 0.33 0.41 0.63
RAN01B 0.24 0.15 0.33
RAN02A 0.44 0.28 0.20
RAN02B 0.31 0.29 0.24
RAN03A 0.01 0.37 0.45
RAN03B -0.14 0.58 0.53
RAN04A -0.32 0.25 0.17
RAN05A 0.13 0.32 0.51
RAN05B 0.14 0.48 0.67
RAN06A -0.10 0.48 0.47
RAN07A -0.02 0.08 0.02
RAN07B -0.15 0.09 0.12
RAN08A 0.27 0.39 0.70
RAN08B 0.17 0.48 0.70
RAN09A -0.27 0.24 0.17
RAN09B -0.30 0.29 0.24
RAN10A 0.25 0.36 0.67
RAN11B 1.00 -0.04 0.30
RAN12A -0.04 1.00 0.71
RAN12B 0.30 0.71 1.00
diag(cor_matrix) #La función ¨diag¨ solamente extrae un el vector diagonal de una matriz cuadrada, y en este caso estamos sacando la correlación entre la misma serie. y por eso deberíamos obtener el valor de 1 en toda la diagonal ya que estamos comparando lo mismo con lo mismo.
RAN01A RAN01B RAN02A RAN02B RAN03A RAN03B RAN04A RAN05A RAN05B RAN06A RAN07A RAN07B RAN08A RAN08B RAN09A RAN09B RAN10A RAN11B
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
RAN12A RAN12B
1 1
Hay funciones en el dplR que pueden generar resultados muy similares a los de COFECHA, el cual es un programa que evalúa la sincronicidad de todas las series entre si. para hacer esta evaluación se dividen las series en segmentos correlativos y traslapados, para evaluar las correlaciones de todas las muestras por segmentos y sin tendencias que pueden sesgar los coeficientes de correlación.
A continuación se generaran unos ejemplos usando la función del dplR llamada corr.rwl.seg siempre es bueno poder utilizar los ejemplos que el mismo paquete provee. los cuales pueden acceder utilizando el siguiente comando, ?corr.rwl.seg en donde pueden leer todo sobre la función, y en la parte inferior de la pagina de ayuda verán algunos ejemplos que se pueden ejecutar en la consola.
En los siguientes ejemplos podrás ver que se ejecuta la misma función pero con diferentes argumentos. (los argumentos son los detalles que definirán como se evalúa y ejecuta la función). por ejemplo acá podemos aclarar el largo del segmento que se utiliza para comparar cada serie con el resto de las series seg.length=10, como también podemos indicarle que se remueva toda la autocorrelación que puede existir en las series prewhiten = T, generar un gráfico representando los resultados make.plot = T y delimitar el limite de confianza para su evaluación estadística pcrit = 0.1
nota: los colores rojos en las figuras indican los segmentos que presentan correlaciones por debajo del limite de confianza que se exije.
COFECHA_p_0.1 <- corr.rwl.seg(rwl = Dat_ringwidth, seg.length = 10, prewhiten = T, biweight = T, make.plot = T, pcrit = 0.1)
COFECHA_p_0.1 <- corr.rwl.seg(rwl = Dat_ringwidth, seg.length = 20, prewhiten = T, biweight = T, make.plot = T, pcrit = 0.1)
COFECHA_p_0.05 <- corr.rwl.seg(rwl = Dat_ringwidth, seg.length = 30, prewhiten = T, biweight = T, make.plot = T, pcrit = 0.05)
COFECHA_p_0.1 <- corr.rwl.seg(rwl = Dat_ringwidth, seg.length = 30, prewhiten = T, biweight = T, make.plot = T, pcrit = 0.1)
De todas las opciones que se pueden ejecutar, qué puedes concluir de este tipo de análisis? -que argumentos son importantes? -de que nos sirve este tipo de análisis?
COFECHA_p_0.1
$spearman.rho
1900.1929 1915.1944 1930.1959 1945.1974 1960.1989 1975.2004
RAN01A NA NA NA NA 0.5052280 0.18620690
RAN01B NA NA 0.2556174 0.2142380 0.5924360 0.17775306
RAN02A NA NA NA NA 0.4896552 0.65250278
RAN02B NA NA NA 0.4327030 0.4371524 0.38242492
RAN03A NA NA NA 0.4967742 0.3147942 0.07897664
RAN03B NA NA NA 0.5537264 0.4994438 0.36907675
RAN04A NA NA 0.4251390 0.4909900 0.6275862 0.66451613
RAN05A NA NA 0.1390434 0.2262514 0.4037820 0.39532814
RAN05B NA NA NA 0.3619577 0.5688543 0.53058954
RAN06A NA NA 0.4233593 0.3971079 0.4518354 0.24983315
RAN07A NA NA 0.3294772 0.2840934 0.4371524 0.46963293
RAN07B NA NA 0.5715239 0.6907675 0.3810901 0.26318131
RAN08A NA NA 0.3023359 0.4745273 0.4665184 0.27919911
RAN08B NA 0.5639600 0.3948832 0.3971079 0.6235818 0.52880979
RAN09A NA NA 0.5101224 0.4438265 0.3486096 0.38553949
RAN09B NA NA 0.6200222 0.6796440 0.5595106 0.57196885
RAN10A NA 0.1350389 0.3005562 0.3241379 0.4527253 0.44115684
RAN11B NA NA NA NA 0.3788654 0.44382647
RAN12A NA NA NA 0.3543938 0.3971079 0.49410456
RAN12B NA NA NA 0.5648498 0.4936596 0.37931034
$p.val
1900.1929 1915.1944 1930.1959 1945.1974 1960.1989 1975.2004
RAN01A NA NA NA NA 0.0024487314 1.615356e-01
RAN01B NA NA 0.0861259490 1.272391e-01 0.0003551341 1.729172e-01
RAN02A NA NA NA NA 0.0032976120 6.713042e-05
RAN02B NA NA NA 8.863031e-03 0.0082473992 1.890727e-02
RAN03A NA NA NA 2.882618e-03 0.0452618173 3.386263e-01
RAN03B NA NA NA 8.888079e-04 0.0027389902 2.274730e-02
RAN04A NA NA 0.0099977278 3.216123e-03 0.0001395279 4.608887e-05
RAN05A NA NA 0.2310322792 1.141534e-01 0.0138720501 1.571371e-02
RAN05B NA NA NA 2.503892e-02 0.0006288791 1.466418e-03
RAN06A NA NA 0.0102816158 1.531020e-02 0.0064640336 9.121239e-02
RAN07A NA NA 0.0379430229 6.402816e-02 0.0082473992 4.749131e-03
RAN07B NA NA 0.0005906781 1.913348e-05 0.0192657779 7.978581e-02
RAN08A NA NA 0.0522860139 4.351818e-03 0.0050176884 6.748935e-02
RAN08B NA 0.0007045424 0.0158159438 1.531020e-02 0.0001560500 1.521954e-03
RAN09A NA NA 0.0022241898 7.391411e-03 0.0298326326 1.809184e-02
RAN09B NA NA 0.0001721577 2.804340e-05 0.0007800962 5.845124e-04
RAN10A NA 0.2375856870 0.0533536507 4.048901e-02 0.0063673391 7.724380e-03
RAN11B NA NA NA NA 0.0198755491 7.391411e-03
RAN12A NA NA NA 2.767294e-02 0.0153102020 3.032660e-03
RAN12B NA NA NA 6.902227e-04 0.0030583065 1.975236e-02
$overall
rho p-val
RAN01A 0.4687628 1.571104e-04
RAN01B 0.2615542 1.091805e-02
RAN02A 0.4532840 7.664083e-04
RAN02B 0.4802881 1.240420e-05
RAN03A 0.4500615 2.722753e-05
RAN03B 0.5190476 7.485563e-06
RAN04A 0.5716716 1.240727e-08
RAN05A 0.3040876 2.980242e-03
RAN05B 0.4712274 2.077248e-05
RAN06A 0.4509711 1.461786e-05
RAN07A 0.3609079 5.009908e-04
RAN07B 0.4990158 6.170749e-07
RAN08A 0.3769783 2.426000e-04
RAN08B 0.5247908 6.316623e-08
RAN09A 0.4878802 6.511242e-07
RAN09B 0.6065775 0.000000e+00
RAN10A 0.3606800 2.242544e-04
RAN11B 0.3299825 5.507254e-03
RAN12A 0.4130697 2.124918e-04
RAN12B 0.5122037 3.975819e-06
$avg.seg.rho
1900.1929 1915.1944 1930.1959 1945.1974 1960.1989 1975.2004
NaN 0.3494994 0.3883709 0.4345351 0.4714794 0.3971969
$flags
RAN01A RAN01B RAN03A RAN05A RAN10A
"1975.2004" "1945.1974, 1975.2004" "1975.2004" "1930.1959, 1945.1974" "1915.1944"
$bins
[,1] [,2]
[1,] 1900 1929
[2,] 1915 1944
[3,] 1930 1959
[4,] 1945 1974
[5,] 1960 1989
[6,] 1975 2004
$rwi
RAN01A RAN01B RAN02A RAN02B RAN03A RAN03B RAN04A RAN05A RAN05B RAN06A RAN07A RAN07B
1865 NA NA NA NA NA NA NA NA NA NA NA NA
1866 NA NA NA NA NA NA NA NA NA NA NA NA
1867 NA NA NA NA NA NA NA NA NA NA NA NA
1868 NA NA NA NA NA NA NA NA NA NA NA NA
1869 NA NA NA NA NA NA NA NA NA NA NA NA
1870 NA NA NA NA NA NA NA NA NA NA NA NA
1871 NA NA NA NA NA NA NA NA NA NA NA NA
1872 NA NA NA NA NA NA NA NA NA NA NA NA
1873 NA NA NA NA NA NA NA NA NA NA NA NA
1874 NA NA NA NA NA NA NA NA NA NA NA NA
1875 NA NA NA NA NA NA NA NA NA NA NA NA
1876 NA NA NA NA NA NA NA NA NA NA NA NA
1877 NA NA NA NA NA NA NA NA NA NA NA NA
1878 NA NA NA NA NA NA NA NA NA NA NA NA
1879 NA NA NA NA NA NA NA NA NA NA NA NA
1880 NA NA NA NA NA NA NA NA NA NA NA NA
1881 NA NA NA NA NA NA NA NA NA NA NA NA
1882 NA NA NA NA NA NA NA NA NA NA NA NA
1883 NA NA NA NA NA NA NA NA NA NA NA NA
1884 NA NA NA NA NA NA NA NA NA NA NA NA
1885 NA NA NA NA NA NA NA NA NA NA NA NA
1886 NA NA NA NA NA NA NA NA NA NA NA NA
1887 NA NA NA NA NA NA NA NA NA NA NA NA
1888 NA NA NA NA NA NA NA NA NA NA NA NA
1889 NA NA NA NA NA NA NA NA NA NA NA NA
1890 NA NA NA NA NA NA NA NA NA NA NA NA
1891 NA NA NA NA NA NA NA NA NA NA NA NA
1892 NA NA NA NA NA NA NA NA NA NA NA NA
1893 NA NA NA NA NA NA NA NA NA NA NA NA
1894 NA NA NA NA NA NA NA NA NA NA NA NA
1895 NA NA NA NA NA NA NA NA NA NA NA NA
1896 NA NA NA NA NA NA NA NA NA NA NA NA
1897 NA NA NA NA NA NA NA NA NA NA NA NA
1898 NA NA NA NA NA NA NA NA NA NA NA NA
1899 NA NA NA NA NA NA NA NA NA NA NA NA
1900 NA NA NA NA NA NA NA NA NA NA NA NA
1901 NA NA NA NA NA NA NA NA NA NA NA NA
1902 NA NA NA NA NA NA NA NA NA NA NA NA
1903 NA NA NA NA NA NA NA NA NA NA NA NA
1904 NA NA NA NA NA NA NA NA NA NA NA NA
1905 NA NA NA NA NA NA NA NA NA NA NA NA
1906 NA NA NA NA NA NA NA NA NA NA NA NA
1907 NA NA NA NA NA NA NA NA NA NA NA NA
1908 NA NA NA NA NA NA NA NA NA NA NA NA
1909 NA NA NA NA NA NA NA NA NA NA NA NA
1910 NA NA NA NA NA NA NA NA NA NA NA NA
1911 NA NA NA NA NA NA NA NA NA NA NA NA
1912 NA NA NA NA NA NA NA NA NA NA NA NA
1913 NA NA NA NA NA NA NA NA NA NA NA NA
1914 NA NA NA NA NA NA NA NA NA NA NA NA
RAN08A RAN08B RAN09A RAN09B RAN10A RAN11B RAN12A RAN12B
1865 NA NA NA NA NA NA NA NA
1866 NA NA NA NA NA NA NA NA
1867 NA NA NA NA NA NA NA NA
1868 NA NA NA NA NA NA NA NA
1869 NA NA NA NA NA NA NA NA
1870 NA NA NA NA NA NA NA NA
1871 NA NA NA NA NA NA NA NA
1872 NA NA NA NA 0.94651776 NA NA NA
1873 NA NA NA NA 0.90397343 NA NA NA
1874 NA NA NA NA 0.57242195 NA NA NA
1875 NA NA NA NA 1.12976698 NA NA NA
1876 NA NA NA NA 0.76456242 NA NA NA
1877 NA NA NA NA 1.14139152 NA NA NA
1878 NA NA NA NA 0.72652422 NA NA NA
1879 NA NA NA NA 1.05778896 NA NA NA
1880 NA NA NA NA 1.01903887 NA NA NA
1881 NA NA NA NA 1.33742461 NA NA NA
1882 NA NA NA NA 0.75486804 NA NA NA
1883 NA NA NA NA 1.10173068 NA NA NA
1884 NA NA NA NA 0.83532752 NA NA NA
1885 NA NA NA NA 0.84445707 NA NA NA
1886 NA NA NA NA 1.07463069 NA NA NA
1887 NA NA NA NA 0.80557260 NA NA NA
1888 NA NA NA NA 1.08410699 NA NA NA
1889 NA NA NA NA 0.85647449 NA NA NA
1890 NA NA NA NA 1.14562403 NA NA NA
1891 NA NA NA NA 0.60590956 NA NA NA
1892 NA NA NA NA 1.13200724 NA NA NA
1893 NA NA NA NA 0.85708224 NA NA NA
1894 NA NA NA NA 0.84641965 NA NA NA
1895 NA NA NA NA 1.19025571 NA NA NA
1896 NA NA NA NA 0.96574328 NA NA NA
1897 NA NA NA NA 1.41057017 NA NA NA
1898 NA NA NA NA 1.07557579 NA NA NA
1899 NA NA NA NA 0.76823276 NA NA NA
1900 NA NA NA NA 1.27857764 NA NA NA
1901 NA NA NA NA 0.75577425 NA NA NA
1902 NA NA NA NA 1.31200133 NA NA NA
1903 NA NA NA NA 0.71761940 NA NA NA
1904 NA NA NA NA 0.85518087 NA NA NA
1905 NA NA NA NA 0.88270042 NA NA NA
1906 NA NA NA NA 0.92232665 NA NA NA
1907 NA NA NA NA 0.98086993 NA NA NA
1908 NA NA NA NA 1.19019025 NA NA NA
1909 NA NA NA NA 0.83568290 NA NA NA
1910 NA NA NA NA 0.89317761 NA NA NA
1911 NA NA NA NA 1.13548355 NA NA NA
1912 NA NA NA NA 0.89662268 NA NA NA
1913 NA NA NA NA 0.95042806 NA NA NA
1914 NA NA NA NA 0.92426289 NA NA NA
[ reached getOption("max.print") -- omitted 92 rows ]
$seg.lag
[1] 15
$seg.length
[1] 30
$pcrit
[1] 0.1
$label.cex
[1] 1
attr(,"class")
[1] "list" "crs"
Como hemos podido observar hay segmentos que no correlacionan con el conjunto de series. Por lo tanto podemos buscar los segmentos que no correlacionan con el resto de series utilizando el objeto resultado de la función corr.rwl.seg el cual se llama COFECHA_p_0.1y podemos indagar adentro de ese objeto utilizando COFECHA_p_0.1$flags.
COFECHA_p_0.1$flags
RAN01A RAN01B RAN03A RAN05A RAN10A
"1975.2004" "1945.1974, 1975.2004" "1975.2004" "1930.1959, 1945.1974" "1915.1944"
#se puede especificar la muestra que te interesa corregir
COFECHA_p_0.1$flags["RAN10A"]
RAN10A
"1915.1944"
COFECHA_p_0.1$flags["RAN05A"]
RAN05A
"1930.1959, 1945.1974"
# Graficar la muestra problemática con otras muestras que se ven bien
# En dplR los nombres de las columnas son las muestras y el nombre de las filas son los años
# Acá podemos extraer el año de inicio (el año mas antiguo del archivo que tenemos)
Año_inicio <- as.numeric(rownames(Dat_ringwidth))[1]
time_series_convert <- ts(Dat_ringwidth[,c("RAN08B","RAN10A")],start = Año_inicio, frequency =1 )
plot(time_series_convert)
# Selección de un solo segmento de las series
segmento_Problemas <- window(time_series_convert, start=1915, end=1944)
segmento_SIN_Problemas <- window(time_series_convert,1970,1990)
# visualizamos el segmento con problemas con respecto a una seria sin problemas
plot(segmento_Problemas[,"RAN08B"])
lines(segmento_Problemas[,"RAN10A"], lwd=2)
# Tambien se puede identificar la columna por su numero
plot(segmento_SIN_Problemas[,1]) #"RAN08B"
lines(segmento_SIN_Problemas[,2], lwd=2) #"RAN10A"
# Ajustamos los limites del Eje "Y"
plot(segmento_SIN_Problemas[,1], ylim= range(segmento_SIN_Problemas))
lines(segmento_SIN_Problemas[,2], lwd=2)
Primero podremos fabricar una serie de tiempo que se apegue a la teoría de el modelo linear agregado “Linear aggregated model” en donde se suman distintas fuentes de variabilidad como el error estocástico, random, Disturbios endógenos y exógenos, Variabilidad climática tendencias relacionadas con la edad.
https://drive.google.com/uc?export=download&id=1yk4RO9h6o9Xd-dbryFErCsQ302xW-9mf
# Acá vamos a generar y crear diferentes vectores de de 200 años, con diferentes magnitudes, tendencias y variabilidad.
Edad <- 1*(1/1.04)^c(1:200) # formula de una curva exponencial negativa
clima <- rnorm(200,mean= 1, sd = 0.3) # se generara un vector con variabilidad aleatoria con media de uno y desviación de 0.3
disturbios <- c(rep(x = 0, 100), seq(from = 0.8,to = 0.001, length.out = 30),rep(0,70)) #disturbio agregado en el año 101 con efecto de 30 años
Error <- rnorm(200, mean = 0.3,sd = 0.1 )
# luego sumamos cada uno de estos factores.
Agregado <- Error+disturbios+clima+Edad
El layout es útil para hacer plots en la misma ventana esta basado en una matriz en donde los números correlativos 1…x es el orden y posición en que hacemos el plot
layout(matrix(1:2, nrow = 2, ncol = 1),widths = 3, heights = c(2,2), respect = T )
par(mar=c(2,4,3,2))
plot.ts(Error, ylim=c(0,2), ylab="variables")
lines(disturbios)
lines(clima)
lines(Edad)
plot.ts(Agregado, main="variabilidad agregada
Error+clima+edad+disturbio")
NA
NA
#Otro ejmplo de uso del layout es:
layout(matrix(c(1,0,
2,5,
3,5,
4,0), nrow = 4, ncol = 2, byrow = T),widths = c(2,3), heights = c(1,1), respect = T )
layout.show(5) # esto solo enseña cual es la configuración de los 5 plots.
par(mar=c(1,4,1,1))
plot.ts(Error, ylim=c(0,2))
plot.ts(disturbios, ylim=c(0,2))
plot.ts(clima, col="darkblue", ylim=c(0,2))
plot.ts(Edad,col="darkgreen", lwd=3, ylim=c(0,2))
par(mar=c(2,4,3,2))
plot.ts(Agregado, main="variabilidad agregada
Error+clima+edad+disturbio")
Primero revisamos la correlación que excite en nuestra seria agregada versusla variabilidad del clima
plot(Agregado,clima) # señal del "agregado" versus el clima
cor(Agregado,clima) # señal del "agregado" versus el clima
[1] 0.74549
Como podemos observar esta relación contiene variabilidad ajena al clima. Recordemos que el agregado contiene error random, un disturbio, y el efecto de la edad, por lo que la correlación no es perfecta y se puede mejorar utilizando algunas tecnicas que a continuación se realizaran.
En los siguientes pasos vamos a experimentar algunas técnicas para remover el efecto de la edad, y el efecto del disturbio utilizando herramientas de estadística y de moldeamiento linear.
Las funciones que se pueden utilizar para esta remoción de tendencias son: detrend.series la que funciona para remover la tendencia a series en especifico y detrend para remover la tendencia simultáneamente a muchas series.
detrend.series(Agregado, method = "ModNegExp") # en este caso utilizamos la función con el argumento del método el cual se utiliza para el ajuste linear. y se selecciona el Modo exponencial negativo para que el programa ajuste una función exponencial negativa a la serie del agregado y se puedan estimar los residuales de este ajuste, para luego estandarizar con una media de uno.
[1] 0.9650803 0.8184300 1.0846493 1.0619995 1.0444616 0.8324715 1.2450632 1.1053481 1.0702778 1.4199945 0.5346874 1.1503803
[13] 0.6851386 1.0592525 0.7822118 0.8668046 0.9949907 1.1587660 1.0801909 0.9179312 1.1670731 0.7420182 1.1168609 0.8877511
[25] 1.2351246 1.3347559 0.9852986 1.1606264 0.9283257 0.8584937 0.9599213 0.8933626 0.8653414 0.7310046 1.4613913 0.9726785
[37] 1.0007010 1.2974423 0.8350237 1.1114405 0.8953357 0.6652086 0.9612884 1.0097921 1.1852772 0.8755057 0.8029161 1.1157914
[49] 1.2075357 0.7839757 1.2618393 0.9550957 0.8888642 0.9526167 0.8888634 0.7267770 0.9663557 1.2523123 0.9671283 0.7379651
[61] 1.3765838 0.8190412 0.9334532 0.8045313 0.9836027 0.7787473 1.2291341 0.8064721 0.6598078 0.8818028 0.5702089 0.8192713
[73] 1.0512476 0.8112225 0.7217231 0.8669991 1.3365290 1.1716990 0.8178966 1.0435956 1.0836355 1.1213571 0.4353135 0.6899652
[85] 1.0368307 1.2683547 0.8690056 1.0711653 0.7988278 0.7829049 1.5457736 1.1485260 0.6663621 1.3583201 1.3776288 1.0626969
[97] 0.8203470 0.7249935 0.8085477 0.8790158 1.7493363 1.6004425 1.2169673 1.4227351 1.4951656 1.0129786 1.1280204 1.4376073
[109] 1.3340330 1.1114918 1.1866664 1.2821141 1.3125411 1.1900140 1.6816308 1.2968492 1.1359497 1.5796055 1.0270297 0.8396560
[121] 1.0356572 1.1134105 0.8112354 1.0756109 1.1467840 1.2917676 0.3730837 0.9367120 0.9921733 0.7428575 1.1082001 0.8121328
[133] 1.1141850 0.8697468 1.1001758 0.7342763 0.7110900 0.8938074 1.1699020 0.5406334 1.3681548 0.9595076 1.3282326 1.2885514
[145] 0.9272305 1.2959988 0.8164180 0.9575205 0.7378651 0.3564974 1.1471446 1.0902738 0.9755660 0.7920068 1.0584035 0.8766122
[157] 1.0131137 1.2702639 0.9066680 0.5979569 0.8394184 0.9794224 0.9948159 1.4979649 1.0764970 0.7022902 0.8531726 0.9653034
[169] 1.0403364 1.1521023 0.5516231 0.5773804 1.2813543 1.1041211 0.8845265 0.7016756 0.8467478 1.0080475 1.0792277 1.1040538
[181] 1.0783939 0.9403540 1.0975040 0.8186996 0.6008448 0.8817551 0.8687051 1.1913259 1.0102961 0.5578136 0.6622438 1.1957642
[193] 0.4872399 1.2236182 1.1028978 1.1756955 0.6753024 1.2224177 1.0545885 0.9739448
detrend.series(Agregado, method = "Spline", nyrs = 20) # en ese caso se ajusta una curva mas flexible llamada Qubic spline con 40 años de flexibilidad.
[1] 1.0083484 0.8383409 1.0909502 1.0516546 1.0217467 0.8074038 1.2019905 1.0677889 1.0403619 1.3960046 0.5335009 1.1645922
[13] 0.7020468 1.0936868 0.8096645 0.8945830 1.0193398 1.1757205 1.0851180 0.9132834 1.1504114 0.7250016 1.0819390 0.8538806
[25] 1.1828352 1.2786513 0.9491423 1.1285664 0.9130119 0.8535198 0.9621707 0.8994532 0.8715875 0.7340471 1.4610670 0.9701226
[37] 0.9984259 1.2988852 0.8412260 1.1279672 0.9149203 0.6827203 0.9871759 1.0347646 1.2106130 0.8912217 0.8143676 1.1278335
[49] 1.2186387 0.7919947 1.2787219 0.9729635 0.9109253 0.9813034 0.9187322 0.7519972 0.9990548 1.2939721 1.0011481 0.7673458
[61] 1.4417742 0.8672829 1.0017796 0.8761781 1.0874716 0.8742376 1.4008166 0.9325380 0.7715026 1.0363399 0.6687225 0.9515862
[73] 1.2025821 0.9111530 0.7944731 0.9353556 1.4176164 1.2294538 0.8543285 1.0900364 1.1354926 1.1805047 0.4596353 0.7267519
[85] 1.0836327 1.3116343 0.8883873 1.0814600 0.7956946 0.7686013 1.4964373 1.0996407 0.6324466 1.2791223 1.2888317 0.9876920
[97] 0.7547748 0.6557908 0.7134241 0.7523774 1.4513747 1.2932192 0.9642623 1.1124435 1.1601402 0.7831503 0.8704230 1.1080547
[109] 1.0279519 0.8566071 0.9145554 0.9884355 1.0138135 0.9237259 1.3178738 1.0325403 0.9244408 1.3208726 0.8863432 0.7490504
[121] 0.9540686 1.0571608 0.7920929 1.0772400 1.1760606 1.3549183 0.3992325 1.0168808 1.0866906 0.8171757 1.2198165 0.8924890
[133] 1.2204424 0.9485644 1.1926651 0.7891692 0.7544592 0.9320694 1.1961091 0.5418263 1.3461672 0.9316804 1.2820760 1.2474794
[145] 0.9081385 1.2925492 0.8326438 0.9987761 0.7843301 0.3834747 1.2387331 1.1761347 1.0487416 0.8470563 1.1248297 0.9257154
[157] 1.0637622 1.3281523 0.9454278 0.6214666 0.8676040 1.0055995 1.0160925 1.5282787 1.1034090 0.7260488 0.8904740 1.0168261
[169] 1.1054856 1.2341086 0.5944016 0.6230285 1.3789600 1.1837313 0.9444294 0.7452339 0.8931114 1.0557882 1.1250489 1.1505695
[181] 1.1294502 0.9946718 1.1765241 0.8908188 0.6625853 0.9818039 0.9731984 1.3392813 1.1374542 0.6266065 0.7376100 1.3126158
[193] 0.5250383 1.2907400 1.1392675 1.1915809 0.6727271 1.1978361 1.0177294 0.9265167
# Acá repetimos lo mismo pero le asignamos un nombre a cada objeto nuevo
det_exp <- detrend.series(Agregado, method = "ModNegExp")
det_sp <- detrend.series(Agregado, method = "Spline", nyrs = 20)
layout(matrix(1:2, nrow = 2, ncol = 1),widths = 3, heights = c(2,2), respect = T )
par(mar=c(2,4,3,2))
plot.ts(det_exp)
lines(det_sp, col="red")
abline(h=1)
plot(det_exp,clima)
points(det_sp,clima, col="red", pch=19)
text(x = 0.6, y = 1.5,labels = round(cor(det_exp,clima),2))
text(x = 0.6, y = 1,labels = round(cor(det_sp,clima),2), col= "red")
cor(det_sp,Error)
[1] 0.379897
require(graphics)
m <- decompose(co2)
plot(m)
Luego de remover las tendencias ajenas al clima procedemos a generar la cronología maestra utilizando todas las series correctamente fechadas. para generar la cronología media de anchos de anillos podemos utilizar la función llamada chron Esta función construye una cronología del valor medio de todas las serias previamente estandarizadas. Típicamente a partir de un data.frame de anchos de anillo sin tendencia según lo producido por detrend.
# Recordamos que tipo de datos tenemos
rwl.report(Dat_ringwidth)
Number of dated series: 20
Number of measurements: 1618
Avg series length: 80.9
Range: 142
Span: 1865 - 2006
Mean (Std dev) series intercorrelation: 0.4451521 (0.08958034)
Mean (Std dev) AR1: 0.7406 (0.08279931)
-------------
Years with absent rings listed by series
None
-------------
Years with internal NA values listed by series
None
# Podemos explotar las estadísticas descriptivas de la base de datos
rwl.stats(Dat_ringwidth)
# Luego pasamos a remover las tendencias de la edad y de posibles disturbios
Dat_ringwidth_index <- detrend(Dat_ringwidth,method = "Spline", nyrs = 50, make.plot = T)
NA
Statistics <- rwi.stats.running (Dat_ringwidth_index, window.length = 20)
Dat_ringwidth_cronologia <- chron(Dat_ringwidth_index, biweight = T,prewhiten = T)
plot(Dat_ringwidth_cronologia)